#pragma once
extern "C" {
#include <f2c.h>
#include <clapack.h>
}
#include "IComplexMatrix.h"

namespace LatoolNet {
	using namespace System;

	ref class ComplexHermitianBandMatrix : public IComplexMatrix {
	private:
		doublecomplex * ap;
		long * ipiv;
		int m_ku;
		int m_kl;
		int m_rownum;
		int m_colnum;
		bool m_isFactorized;

		ComplexHermitianBandMatrix(ComplexHermitianBandMatrix ^ orig) {
			m_rownum = orig->m_rownum;
			m_colnum = orig->m_colnum;
			
			m_ku = orig->m_ku;
			m_kl = orig->m_kl;


			int size = (m_ku + 1) * m_rownum;
			ap = new doublecomplex[size];

			for (int i = 0; i < size; i++) {
				ap[i] = orig->ap[i];
			}
			ipiv = new long[m_rownum];

			for (int i = 0; i < m_rownum; i++) {
				ipiv[i] = orig->ipiv[i];
			}

			m_isFactorized = orig->m_isFactorized;

		};
	public:
		ComplexHermitianBandMatrix(int rownum, int colnum, int subnum, int supernum) {
			if (rownum != colnum || subnum != supernum) {
				throw gcnew LatoolException("Number of rows and columns must be equal for symmetric matrix.");
			}
			m_rownum = rownum;
			m_colnum = colnum;

			m_ku = supernum;
			m_kl = subnum;

			int size = (supernum + 1) * m_rownum;
			ap = new doublecomplex[size];

			for (int i = 0; i < size; i++) {
				ap[i].r = 0.0;
				ap[i].i = 0.0;
			}

      ipiv = new long[rownum];

			m_isFactorized = false;
		};

		~ComplexHermitianBandMatrix(){
			delete[] ap;
		};

		virtual IComplexMatrix ^ Clone() {		
			return (IComplexMatrix^) gcnew ComplexHermitianBandMatrix(this);
		};

		virtual property int RowNum {
			int get() {
				return m_rownum;
			}
		};
		virtual property int ColNum {
			int get() {
				return m_colnum;
			}
		};
		virtual property Complex default[int, int] {
			Complex get(int i, int j) {
				if (i > j) {
					int temp = i;
					i = j;
					j = temp;
				} 
				if (i == j || j == i + 1) {
					return Complex(ap[(m_ku + 1) * j + m_ku + i - j].r, ap[(m_ku + 1) * j + m_ku + i - j].i);
				} else {
					return Complex(0.0, 0.0);
				}
			}
			void set(int i, int j, Complex value) {
				if (i > j) {
					int temp = i;
					i = j;
					j = temp;
				} 
				if (i == j || j == i + 1) {
					int idx = (m_ku + 1) * j + m_ku + i - j; 
					ap[idx].r = value.Real;
					ap[idx].i = value.Imag;
					m_isFactorized = false;
				}
			}
		};
		virtual property MatrixType Type {
			MatrixType get() { return MatrixType::ComplexHermitianBand; }
		};
		virtual property bool IsFactorized {
			bool get() { return m_isFactorized; }
		}

		virtual void Invert() {
			throw gcnew LatoolException("Inversion of banded matrix is not supported by LAPACK.");
		};

		virtual void Factorize() {

			char uplo = 'U';
			long n = m_rownum;
			long kd = m_ku;
			long ldab = kd + 1;
			long info;

			zpbtrf_(&uplo, &n, &kd, ap, &ldab, &info);

			if (info == 0) {
				m_isFactorized = true;
				return;
			}

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			} else {
				throw gcnew LatoolException("This matrix is not positive definite. Symmetric tridiagonal indefinite matrix is not supported by LAPACK. Please use Matrix::DoubleSymmetrix."); 
			}

		};

		virtual void Solve(IMatrix ^ b) {

			ComplexGeneralMatrix ^ gb = (ComplexGeneralMatrix^) b;

			char uplo = 'U';
			long n = m_rownum;
			long kd = m_ku;
			long nrhs = b->ColNum;
			long ldab = kd + 1;
			long ldb = Math::Max(1, n);
			long info;

			zpbtrs_(&uplo, &n, &kd, &nrhs, ap, &ldab, gb->dat, &ldb, &info);					

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}
		};

		virtual void SolveWithSimpleDriver(IMatrix ^ b) {

			ComplexGeneralMatrix ^ gb = (ComplexGeneralMatrix^) b;

			char uplo = 'U';
			long n = m_rownum;
			long kd = m_ku;
			long nrhs = b->ColNum;
			long ldab = kd + 1;
			long ldb = Math::Max(1, n);
			long info;

			zpbsv_(&uplo, &n, &kd, &nrhs, ap, &ldab, gb->dat, &ldb, &info);
			
			if (info == 0) {
				m_isFactorized = true;
				return;
			}

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			} else {
				throw gcnew LatoolException("This matrix is not positive definite. Hermitian band indefinite matrix is not supported by LAPACK. Please use Matrix::ComplexHermitian."); 
			}
		};


	};
}